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Interacting systems with K driven particle species on a open chain or chains which are coupled at 
the ends to boundary reservoirs with fixed particle densities are considered. We classify discontinuous 
and continuous phase transitions which are driven by adiabatic change of boundary conditions. We 
build minimal paths along which any given boundary driven phase transition (BDPT) is observed 
and reveal kinetic mechanisms governing these transitions. Combining minimal paths, we can drive 
the system from a stationary state with all positive characteristic speeds to a state with all negative 
characteristic speeds, by means of adiabatic changes of the boundary conditions. We show that along 
such composite paths one generically encounters Z discontinuous and 2(K — Z) continuous BDPTs 
with Z taking values < Z < K depending on the path. As model examples we consider solvable 
exclusion processes with product measure states and K = 1, 2, 3 particle species and a non-solvable 
two-way traffic model. Our findings are confirmed by numerical integration of hydrodynamic limit 
equations and by Monte Carlo simulations. Results extend straightforwardly to a wide class of 
driven diffusive systems with several conserved particle species. 



I. INTRODUCTION 

Systems of interacting particles out of equilibrium with 
more than one particle species find wide range of ap- 
plications in biological, social and physical context 
Q, 0| . Remarkable phenomena occurring in these sys- 
tems are the so-called boundary-driven phase transitions 
(BDPTs), e.g. phase transitions induced uniquely by 
changes of the boundary conditions Q. These transi- 
tions can be both of first and second order, depending on 
whether the order parameter (e.g. the stationary bulk 
density) changing discontinuously or continuously across 
the transition line. It has been demonstrated that in sys- 
tems with one particle species the first order (discontinu- 
ous) BDPTs are governed by shocks dynamics, while the 
second order (continuous) BDPTs are governed by rar- 
efaction waves Q . In case of several particle species cor- 
responding in the hydrodynamic limit to the hyperbolic 
systems of conservation laws, it was noted that the num- 
ber of qualitatively different first-order BDPTs increases 
with the number of particle species Q , a fact which is re- 
lated to a complex structure of shocks in systems of con- 
servation laws @, However, generic properties of 
boundary driven phase transitions in multi-species driven 
particle systems are largely unknown. Also, the mecha- 
nisms governing stationary state selection are known for 
one-species systems [H,[Hj, while for systems with several 
particle species such a knowledge is lacking. 

The aim of the present paper is to show the existence 
of hierarchies of qualitatively different phase transitions 
in driven systems with several particle species, with open 
boundaries, and to provide a first classification of them. 
In particular we show that in systems with K species 
of particles there are generically at least K first order 
(discontinuous) and K second order (continuous) quali- 
tatively different phase transitions. We demonstrate that 
these phase transitions can be observed by driving the 
system from a stationary state with all positive char- 



acteristic velocities to a state with all negative charac- 
teristic velocities, by means of adiabatic changes of the 
densities at the boundaries. The kinetic mechanisms un- 
derlying the occurrence of first order BDPTs is ascribed 
to shock waves interactions (both among themselves and 
with boundaries) while the origin of second order phase 
transitions is shown to be connected with rarefaction 
waves. This leads to several physical implications which 
are confirmed by numerical simulations. In particular, 
we show that the existence of shock waves and the con- 
tinuity of the flux across a first order BDPT allows us to 
predict the location of the stationary densities after their 
discontinuous change. Similarly, conditions of stability of 
rarefaction waves allows us to predict the location of the 
second order phase transition points along a specific path 
in parameter space. Also, it follows that the occurence 
of qualitatively different BDPTs is connected to the ex- 
istence of different classes of shock waves and rarefaction 
waves in systems with several particle species. Continu- 
ous paths in parameter space along which sequences of Z 
first order and 2{K — Z) second order transitions with [Z 
= 0, 1, K) occur, are identified. We show that for each 
value of Z, there are (^) qualitatively different paths, 
each of them leading to a different set of transitions. 

Our general approach is tested on ideal solvable mod- 
els (e.g. which admit product steady states) as well as on 
more realistic (not exactly solvable) ones. For the former 
we derive hydrodynamic equations of motion and use the 
theory of shock and rarefaction waves for system of con- 
servation laws for their investigation. For non-solvable 
models we recourse to Monte-Carlo simulations as prin- 
cipal tool of investigation. To simplify the presentation 
we concentrate in more detail on the case of two particle 
species but the results arc discussed in a manner that the 
extension to the case of an arbitrary number of particle 
species becomes straightforward. 

The plan of the paper is the following. In Section |TT] 
we introduce multi-species particle models and their basic 
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properties. In Sec. IIIII we discuss hierarchies of BDPTs 
in multi-species systems and in Sec. IIVI we character- 
ize the "minimal path" along which any chosen (discon- 
tinuous or continuous) BDPT can be observed, taking 
as working examples the cases of one and two particle 
species {K = 1,2). In Sees IVl and IVII we discuss the ba- 
sic mechanisms governing first and second order BDPTs 
in multi-species systems, respectively. In Sec. I Villi we il- 
lustrate our approach for the case of a more realistic (not 
exactly solvable) model for bidirectional traffic, while in 
Sec. IVIII we discuss special paths for models with K 
particle species along which sequences of BDPTs in the 
parameter space occur. Finally, Sec IIXI serves for conclu- 
sions and for open problems. 



II. MULTI-SPECIES PARTICLE MODELS OUT 
OF EQUILIBRIUM 

We consider a system, consisting of many interacting 
particles, dynamics of which is governed by a Master 
Equation, with the rates {Tec} of transition between 
the states C and C", on a discrete lattice. We assume 
that there are K > 1 particle species conserved indepen- 
dently, and that one can construct a hydrodynamic limit, 
which has the form of a system of conservation laws, of 
the type 



dt dx 



-e— (y B k d -^ 
dx \j~ 3 dx 

1,2,...,K. 



(1) 



Here p k (x,t) denotes the averaged local density of the 
specie fc, jk is the respective flux, Bij(pi, p 2 , Pk) is 
the diffusion matrix, e is an infinitcsimally small positive 
constant. The system is defined on a segment x S [0, 1]. 
At the ends of the segment, Dirichlet boundary condi- 
tions are imposed: 



Pk[0] = Pk,L, Pfc[l] = Pk,R- 



(2) 



We are interested in the stationary state of the system, 
i.e. in a state attained in the infinite time limit t —> oo. 
Note that in the case of perfect match between left and 
right boundary densities Eqs. ([I]),© allow stationary 
space-homogeneous solution pk(x) = pk,L = Pk,R for 
all k. Stability of this solution for all physical values 
of Pi, p2i •••) Pk is guaranteed by the positive definitc- 
ness of the diffusion matrix Bij. For the source particle 
model, this corresponds to an existence of a homogeneous 
stationary state with constant average particle densities 

Pl,P2,—,PK- 

Some remarks are in order. Dirichlet boundary condi- 
tions © in terms of particle system mean that the sys- 
tem is coupled to boundary reservoirs with fixed particles 
densities Pk.L-, Pk.R a t the left and right ends, respectively 
How this can be implemented or read off from boundary 
rates was discussed in llj. The flux jk(pii p2, Pk) is a 



nonlinear function of the particle densities, which implies 
interaction between particles. For a precise definition of 
the genuine flux nonlinearity see e.g. Q. 

A special role is played by the Jacobian matrix of the 
flux (P'i) kw = djk/dp w whose eigenvalues c\ < c 2 < ... < 
ck, called characteristic speeds, are all real and distinct 
jl2| (note that we numerate the characteristic speeds in 
the increasing order). The characteristic speeds for the 
original particle system are velocities with which the lo- 
cal perturbations of a homogeneous state are propagating 
[13( | . The physical region (i.e. the region with physically 
meaningful particle densities pk) splits into domains ac- 
cording to the signs of characteristic velocities. E.g. for 

K = 2 we shall call G |_ the domain in u, v space (u, v are 

particle densities of the two species) where c\(u, v) < 
and c 2 {u,v) > 0. Characteristic speeds are smoothly- 
changing functions of particle densities. Special role is 
played by the subdomains of dimension K — 1 across 
which, one of characteristic speeds changes its sign. E.g. 

the domains G++and G |_ are separated by the subdo- 

main Go+ with c\(u, v) = 0. An example of a decompo- 
sition of the physical region for K = 2 is given in FigJU 
If a stationary state of a particle model has bulk parti- 
cle densities u stat ,v s tat belonging to the G++ region, we 
call it G++ stationary state, etc.. Analogously, if the left 

boundary densities belong to, say, G |_ domain, 

we say that left boundary is of G |_ type etc. The gen- 
eralization to arbitrary number of particle species K > 2 
is obvious. 



III. HIERARCHIES OF CONTINUOUS AND 
DISCONTINUOUS BDPT 

Different types of BDPTs can occur in these mod- 
els. A discontinuous (first order) phase transition of the 
p — th type is a transition between a stationary state 
with ci, Cp_i < 0, c p , Ck > and a stationary state 
where c p has changed its sign. E.g., for K = 2, type 1 
and type 2 discontinuous transitions are transitions be- 
tween G ++ ^ G- + and G_ + ^ G__, respectively. Note 
that there is no direct transition between the G++ and 
G__ state. This is due to the strict hyperbolicity: there 
is no region, where G ++ and G__ touch each other (the 
contrary would mean the existence of a weak hyperbolic 
point [12j). Obviously, in systems with K species p can 
generically take values p = 1,2, ...K, leading to K differ- 
ent types of discontinuous phase transitions. 

The continuous (second order) phase transition of type 
p is a transition between a stationary state with zero 
p — th characteristic velocity, c p — 0, and a station- 
ary state where c p is strictly positive or negative. The 
signs of other characteristic velocities c q , q ^ p do not 
change across this transition. E.g. for K = 2, transitions 

G++ ^ Go+,Go+ G |_ are of type p = 1, while the 

transitions G |_ G_o> G-o ^ G__ are of type p = 2. 

Note that continuous transition of type Gq+ ^ G_o does 
not happen since the respective subregions have no inter- 
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Figure 1: A schematic diagram, showing minimal Phi path 
(bold lines) and Ph2 paths (thin lines) in parameter space 
(Panel (a)) and in physical region (Panel (b)). The Ph2 path 
is continuous both in parameter space and in physical region. 
The Phi path is continuous in parameter space but is discon- 
tinuous in physical region. Panel (a): The axes represent 
sets of boundary densities at the left and at the right bound- 
ary. The part of the Ph2 trajectory marked by dashed line 
indicate the rarefaction- wave governed stationary state, (see 
Sec lVI|) . The bold dotted line in the left upper corner marks 
the part of the Phi path, inside which a first order transition 
occurs, see Sec(V] Panel (b): Schematic evolution of the sta- 
tionary densities along a Phi path (bold broken line) and a 
Ph2 path (thin line). Along the Phi path at least one discon- 
tinuous p-type phase transition is observed. Along the Ph2 
path, a pinning and depinning to a state with zero character- 
istic velocity (two continuous phase transitions of p-type) are 
observed. The pinning and depinning point correspond to the 
ends of the dashed-marked Ph2 segment on Panel (a). 



section due to strict hyperbolicity. In systems with weak 
hyperbolic point, where such an intersection exists, a con- 
tinuous transition Gq+ ^ G_o may become possible (see 
direct Ph2 path in FigQT] (b)). In case of K species, p 
can generically take values p = 1,2,. .K, leading to K 
different types of continuous (second order) transitions. 



IV. MINIMAL PATHS TO OBSERVE BDPT 

Here we describe the simplest path in parameter 
space along which a single first or second order phase 
transition is surely observed, in a system with K 
species of particles. For this we introduce a param- 
eter space of dimension 2K, coordinates of which are 
given by the left and right boundary densities {pl\pr} = 
{Pi,l,P2,l, —,Pk,l\pi,r, —,Pk,r}- Each point of the pa- 
rameter space represents some boundary conditions (|2|). 
Let us also define a physical region of dimension K co- 
ordinates of which are average bulk particle densities 
Pl, pi, ..-Pk- Each stationary state with bulk particle 
densities p sta t = {pf at , p|* a * > •■•> PK at } 1S represented by a 
point in the physical region. The bulk stationary density 
will be our order parameter. 

A path in parameter space is represented by a con- 
tinuous curve r(s) given by the left and right boundary 



densities {pl(s)\pr(s)} along the path, parametrized by 
the running coordinate < s < 1, with s = and s = 1 
corresponding to the initial and final points of the path, 
respectively. For each value of s, we wait until the system 
reaches a stationary state (which we assume to be homo- 
geneous in the bulk) and record the stationary bulk par- 
tide densities p stat (s) = {pf at {s), p s 2 tat (s), p K at (s)}. 
Thus, each path T(s) in parameter space of dimension 
2K is mapped on a path p s tat{s) in physical region of 
dimension K. The mapping T(s) — > p s tat(s) is not in- 
vcrtible, because many different boundary conditions can 
lead to states with equal bulk stationary densities. 

Consider two neighboring domains Gx and Gy in the 
physical region p\,p 2 , ■■■Pk characterized by the following 
signs of the characteristic speeds Ci{p\,p2, ...ppc)'- 

G x ■ ci < ... < c p _i < 0; c K > c K -i > ... > c p > 0, 
Gy ■ c\ < ... < c p < 0; ck > ck-i > ■■■ > c p+ i > 0, 

(3) 

Note that both Gx , Gy have dimension K and that do- 
main Gx has one positive characteristic speed more than 
Gy. We denote with Gxoy a subdomain of dimension 
K — 1 separating domains Gx and Gy, given by: 

Gxoy ■ c\ < ■■■ < Cp-i < 0; c p = 0; 
ck > ck-i > ■•■ > c p+ i > 0. 

Phase transitions occur along paths which connect 
neigbouring domains. In the following we denote by 
p™ % , the sets of the particle densities in the left and 
right boundary reservoirs at the initial point s = of 
a path r(s), i.e. = {p*,l(* = 0)}f =1 , p™ = 

{pk, R (s = 0)}f = i, T(s = 0) = {pT 1 \pT}. Analogously, 
denote p^ nal , p^ nal the respective boundary densities at 
the end of the path T(s = 1) = {p^ nal \p^ nal }. Let us 
choose the initial and the final boundary densities on 
both boundaries from domains Gx and Gy, respec- 
tively, 

pt 6 G x ,pt G G x and p{ mal e G Y , P f ™ al e Gy. 

(4) 

In addition, we shall choose the initial and the final 
boundary densities fully matching, = = Pini, 
and p£ ma ' = p-^ nal = ppjj^. Since the full match allows 
for trivial homogeneous stationary solution of (H]),©, the 
initial and final stationary states along the path belong 
to different G domains: p s tat{s = 0) = pini £ Gx and 
Pstat{s = 1) = pfin € Gy, see discussion after Eq.([2]). 
Note that the condition of the full match can be relaxed 
(for an example see Sec lVIIl|) . and is chosen here for sim- 
plicity of presentation. 

Having chosen the initial and final points of a path 
r(s), let us now define paths T phl (s), T ph2 (s) of type 
Phi and Ph2 as consisting of two consecutive elementary 
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Figure 2: Panel (a): Phase diagram of TASEP with open 
boundaries. LD, HD and MC denote Low Density, High Den- 
sity and Maximal Current phases respectively. Thick(thin) 
lines separating different phases mark phase transitions of the 
first ( of the second) order. Thick and thin lines originating 
in INI point (0.3,0.3) and ending in FIN points show Phi 
and Ph2 paths, respectively. Panel (b): Stationary densi- 
ties along the Phi path(bold line) and along the Ph2 path 
(thin line) shown on Panel (a), versus running variable s, 
parametrizing the paths. A cusp in the upper part of the Phi 
density path u„tat(s) results from a cusp in the Phi path in 
parameter space(Panel a). 



steps: 



r™(«) : {pT\pT} {pT\p f R nal } -> {pi mal K nat } 



final I final ^ 



(5) 



(*) : K"Vif } -» {p f L nal \pT} "> {p f L mal \p f H mal }- 



(6) 



During each elementary step the boundary densities 
at one boundary change adiabatically while the other 
boundary is kept fixed. We shall call an elementary 
step during which the left (the right) boundary densities 
change a step L ( step R), respectively. Thus, a Phi path 
T Phl (s) consists of consecutive steps step R — > step L, 
while in a Ph2 path the steps R and L are interchanged, 
see Fig. [TJ It is important that a domain Gx of depar- 
ture for both paths has more more positive characteristic 
velocities than the domain of arrival Gy, the reason for 
which will be explained in the next section. 

We assume trajectories — > p^ nal (step R) and 
pim — j. pf* nal ( s tep L) to lie entirely inside Gx and Gy, 
and to cross c p = hyperplane only once, see FigQ] (a). 
The main results of this section are presented in the two 
propositions listed below. 

Proposition I. Discontinuous change of the station- 
ary density p s t a t from a point in domain Gx to a point 
in neighboring domain Gy (p—th type of first order tran- 
sition) occurs along any Phi path f5|) connecting these 
domains. 

The path T phl (s) is continuous in the parameter space 
coordinates of which are left and right boundary den- 
sities. We claim that at some point s* whose precise 
location depends on the microscopic transition rates, 



the stationary density p s tat(s) undergoes discontinuous 
jump (a first order transition) from some value X G Gx 
(non necessarily coinciding with initial point p™" 1 ) to a 
value Y € Gy (non necessarily coinciding with end point 
pfmaiy rpkg stationary flux across the transition is con- 
tinuous, jstat(X) — Jstat 

Proposition II: Two continuous (second order) phase 
transitions of type p are observed along any path of type 
Ph2. Those transitions are continuous phase transitions 
Gx — > Gxoy and Gxoy — > Gy. 
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Bidirectional traffic model 

Figure 3: Multi-chain model and two-way traffic model on a 
narrow road. Coupling to boundary reservoirs is indicated by 
boxes marked L (left reservoir) and R (right reservoir). 



Summarizing Propositions I and II, we have that by 
proceeding along Phi (Ph2) path connecting two sta- 
tionary states in neighboring regions, we observe one 
first order (two second order) transitions between these 
states (see FigllJ. Alternatively, one can say that any 
Phi path in parameter space maps onto a discontinuous 
path p s tat(s) in physical region, while any Ph2 path in 
parameter space maps onto a continuous path p s t a t (s) in 
the physical region, a segment of which is pinned to the 
c p = hyperplane, see Fig[TJ Continuous phase tran- 
sitions Gx — > Gxoy and Gxoy — > Gy correspond to 
pinning and depinning points of the Ph2 trajectory. We 
remark that in particle systems without hysteresis the 
above paths are fully invertible, i.e. by proceeding along 
a path r(s) in the opposite direction from the final to 
initial point, one encounters exactly the same mapping 
r(s) -> p s tat(s). 

The physical significance of the Ph2 path resides in 
the fact that it contains favourable boundary setups for 
the formation of rarefaction waves governing the Gxoy- 
type stationary states (Ph2 segment marked in FigfTJa) 
by dashed line) and does not contain boundary setups 
favoring stable shock waves which govern discontinuous 
phase transitions. Consequently, discontinuous shock- 
wave driven transitions are impossible along any Ph2 
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Phi 

path and the mapping T(s) — > p s tat(s) is continuous (see 
also sec lVI[) . Analogously, one can deduce that the map- 

Phl 

ping r(s) — > Pstat(s) along any Phi path must be dis- 
continuous. More details are given in Secs lVlVTl Before 
discussing the kinetic mechanisms underlying the transi- 
tions we illustrate Propositions I and II with some specific 
examples. 

Case K=l. It is instructive to start with the sim- 
plest possible and well known case of one particle species 
K = 1. Let us take the Totally Asymmetric Simple Ex- 
clusion Process, or TASEP [l4|],[l5[ as a representative. 
This process is defined on a chain, each site of which can 
be empty or occupied by one particle. Particles jump 
independently after an exponentially distributed random 
time with mean 1 to a nearest neighbor site on the right, 
provided that the target site is empty (hard core exclu- 
sion rule). We recall that the particle flux, which is a 
number of particles crossing a single bond per unit time, 
as function of the density u for this model has the form 
j(u) = u(l — u) and the physical region of densities is 
< u < 1 due to the hard core exclusion. The character- 



Case K=2. As examples we choose a multi-chain 
model [l6[ restricted to K = 2, and a two-way traffic 
model on a narrow road. The former describes interact- 
ing exclusion processes evolving on parallel chains with 
species hopping in the same direction, while the latter 
describes the exclusion process with species hopping in 
opposite direction, see Figf3] Both models are described 
in continuous limit by equations of type (H},©. 

The multi-chain model is solvable on a ring, its sta- 
tionary current is known analytically together with the 
diffusion matrix B, this giving the possibility of using ei- 
ther hydrodynamic limit equations ([1]) or microscopic ap- 
proach (Monte Carlo simulations). Moreover, it is known 
from previous studies [13], [l6|, that numerical integra- 
tions of the respective discrctized hydrodynamic equa- 
tions ([1]), ©, reproduce very accurately (if not exactly) 
the outcome of the respective Monte- Carlo simulations. 
The other model (two way traffic) is not solvable and 
Monte Carlo simulations remain the only tool of investi- 
gation. 
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Figure 4: Decomposition of the physical domain for a model 
with two species [ multi-chain model with K = 2, 7 = 0.5, see 
text for model definition], according to signs of characteristic 
speeds ci,C2. E.g. Go+ denotes region where ci = 0, ci > 
etc.. Points A,B,C,D are reference points in respective 
G-domains, chosen for illustration of Phi and Ph2 paths in 
FigsEU 

istic speed c = j'(u) = 1 — 2u is positive for u G [0, 1/2), 
is negative for u E (1/2, 1], and it vanishes for u = 1/2, 
thus defining the respective domains G+,G_ and Go- 
The boundary densities arc ul and ur. The well-known 
stationary states of TASEP with open boundaries, the 
Low Density (LD) state, the High Density (HD) state, 
and the Maximal Current (MC) state are readily identi- 
fied with G+, G_ and Go states, respectively. The phase 
diagram of TASEP and the mappings T(s) — > u sta t(s) 
for Phi- and Ph2-paths are shown in FigJ5] Note that 
along the Phi path one encounters a discontinuous phase 
transition from LD to HD state G+ — > G_. Along the 
Ph2 path two continuous phase transitions G_|_ — > Go , 
Go — ^ G_ take place. 
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Figure 5: Location of stationary densities for two-chain model 
(7 = 0.5) for (a) two alternative Phi paths (b) two alter- 
native Ph2 paths from G++ to G ^ domains. Initial and 

end points for all paths corresponds to fully matching left- 
right boundary densities p ini = (0.2,0.3), pf ina i = (0.6,0.4) 
(points A,B in Fig(4]). Path I (filled symbols) is a direct path 
Pini — > pfinai, where the boundary densities during respec- 
tive steps L and R (see Sec lIV|) change by linear interpola- 
tion, e.g. u R (s) - i4T = ("fl" 1 - v?r % )s, etc.. Path II (open 
symbols) goes through intermediate point C with coordinates 
pc = (0.3, 0.5) i.e. each step L,R involves sequence of two lin- 
ear interpolations — ► pc,pc pfinai- Evolution of the 
densities along the paths I and II is indicated by solid and 
dashed arrows, respectively, crosses mark initial, intermedi- 
ate and final points. Data points are taken from numerical 
integration of the respective discretized hydrodynamic equa- 
tions Q, Deviations from the exact c p — line in Panel 
(b) are due to finite-size errors. 

First consider the multi-lane model for K = 2, see the 
top panel of Fig[3J The only allowed move is a hopping 
of a particle to its nearest neighbouring site on its right, 
with the rate r n = 1 — 717/2 where < n < 2 is the num- 
ber of particles on the adjacent chain, neigbouring to the 
departure and to the target site (see Fig[3|). The inter- 
chain interaction parameter 7 varies from 7 = (corre- 



G 



sponding to two uncoupled TASEPs ) to 7 = 1 (maximal 
interaction). Particle hopping obeys the exclusion rule: 
if the target site is already occupied, the move is rejected. 
The hopping between chains is not allowed. We denote 
the average density of particles on the first and second 
chain as u, v, respectively. Due to hardcore exclusion, the 
physical region of densities for this model is the square 
domain < u, v < 1. The model has product stationary 
states, which allows to calculate stationary currents [16[ 
in = u (l — — l v ) an d jv = w (l — — l u )- The 
characteristic velocities are the eigenvalues of the flux 
Jacobian (Dj)ipk = c^tpk, where 



(1 



- 2u)(l 
-7^(1 — 



(1 



-7it(l — u) 
- 2v)(l -711) 



(7) 



Consistently with our notations, a domain in u,d plane 
with C\{u,v) > 0, 02(11, v) > will be denoted by 
G++ (and similarly for other sign combinations). For 

generic interaction 7 < 1, all domains G++,G |_,G 

are present, see Fig|4j As 7 increases, the domain G 

shrinks and for 7 = 1 (maximal interaction) it disappears 
completely. 

According to the Proposition I, described in the Scc lIVl 
starting from the domain G++ and going to the domain 
G |_ using a Phi path ©, we should observe discontin- 
uous transition in stationary density along the path. In 
Fig. 5(a) stationary densities along two alternative Phi 
paths between the same initial point (located in G++) 

and final point (located in G |_) are shown. These paths, 

indicated as AB and ACB in FigJU differ by the elemen- 
tary trajectories pi n i p final are built, see caption of 
Fig 5(a) As expected, the qualitative scenario does not 
depend on the details of the path: in both cases we see a 

discontinuous transition G++ — > G |_. Similarly, along 

the two alternative Ph2 paths © two continuous BDPT 
occur: G ++ — > Go+ and Go+ — > G_ + (sec Fig 5(b)), in 
accordance with the Proposition 2. 
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Figure 6: Panels (a)-(d): Location of stationary 
densities for two-chain model with 7 = 0.5 along 
(Phl)(Phl),(Ph2)(Ph2),(Phl)(Ph2), and (Ph2)(Phl) paths, 
respectively. All trajectories pass through initial (A), middle 
(B) and final (D) points. E.g. Panel (b) shows stationary den- 
sities along two consecutive Ph paths A ~^ phl B,B — > ph2 D. 
Filled (empty) symbols mark stationary densities u s tat, v s tat 
along the first (the second) Ph path. Evolution of the den- 
sities along the path is shown by arrows. Along every Ph2 
piece a pinning to/depinning from c p — line is observed 
(two continuous phase transitions). Along every Phi piece 
a discontinuous phase transition occurs, marked by dotted 
arrows. 



A. Building composite paths 



B. BDPTs for arbitrary number of species 



From the point B in G |_ where all above described 

minimal paths ended, we can continue further building 

Phi G_+ ^ GL_ © or Ph2 G_+ P -f G__ © paths 
to some point D in G__, see Fig 01 Along the new Phi 
(Ph2) path from B to D we shall see another discon- 
tinuous (continuous) transition in the stationary density. 
Combining all possible Phi and Ph2 paths leading from 
G-\ — l — ^ G — l — y G — , one can observe a desired se- 
quence of phase transitions. E.g. choosing only Phi 



G , we see two phase tran- 

Along a path 



^ Phi ~ Phi 

paths G ++ — > G_ + -> 
sitions of the first order, see Fig 6(a) 

Ph2 Phi 

G ++ -» G_ + — > G__ we see two continuous transi- 
tions G++ — > Go+ — > G |_ and a discontinuous transi- 
tion G |_ —¥ G__, etc.. The outcome of all four possible 

choices of composite paths passing from G ++ to G are 

shown in Figs 6(a)|6(dJ 



The generalization of these considerations to systems 
with arbitrary number of species is straightforward. In 
systems with K particle species, one can construct com- 
posite paths leading from a state with all positive char- 
acteristic speeds G+...+ to a state with all negative char- 
acteristic speeds G_..._, similarly to those shown in 
Figs 6(a)||6(d) Every composite path consists of K con- 
secutive minimal Phi or Ph2 paths, described in Scc lIVl 
Obviously, the number of qualitatively different compos- 
ite paths is 2 K . Along any composite path from G+...+ 
state to G_..._ state we shall observe at least Z first or- 
der transitions, and at least 2{K—Z) discontinuous phase 
transitions, where Z is the number of Phi pieces in the 
composite path. We say "at least" , because e.g. the num- 
ber of first order transitions can occasionally be larger 
than Z (in case of complicated shape of G-domains or if 
composite path crosses the same hypcrplane c„ = sev- 
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Figure 7: Average density profiles u(x,t),v(x,t) evolution 
close to the first order phase transition line. Initial state at 
t = is a state matching the left boundary ul = 0.36, vl = 
0.34(the right boundary is ur — 0.6, vr = 0.4). As time goes 
on, a shock wave at the right boundary appears and starts to 
propagate, its position at t = 6400 is shown. After reaching 
the left boundary, it reflects, see Inset, taken at t — 17600. 
Final stationary state has densities u = 0.5568, v — 0.4868. 



cral times), but not smaller than Z, and analogously for 
second order transitions. E.g. already in one-component 



multiple reflections 
reflected state 




space 

Figure 8: Schematic picture of a shock waves interaction 
scenario for two-species system. Dashed line denote multi- 
ple reflections of the — + shock from the boundaries, during 
which the bulk density exponentially converges to its station- 
ary value. 

systems with a double maximum in the current-density 
relation j(u), both G+ and G_ domains consist of two 
separated segments and a Phi path between disjoint G+ 
and G— domains may lead to observation of two discon- 
tinuous transitions. In the next two sections we discuss 



the kinetic mechanisms governing BDPTs along minimal 
Phi and Ph2 paths. 



V. SHOCK WAVE MECHANISM UNDERLYING 
FIRST ORDER BDPT 

Consider the first order phase transitions shown in 

Fig 5(a) from the G ++ to G |_ state. Initial G ++ state 

is a state perfectly matching the left boundary. Note 
that the perfect match with the left boundary is the only 
possible way for a G++ state to be stationary because 
any eventual perturbation at the left boundary will be 
carried away from it since both characteristic velocities 
are positive (see also |f|). As we cross a transition point 
(at a small but finite distance from it) a new — h shock 

wave with densities u\ + , vf + belonging to G [_ appears 

at the right boundary and starts to propagate inside the 
bulk. Note that the new shock does not match the right 
boundary ^ ujj,v^ + ^ Vr) but forms a bound- 

ary layer with it (see FigEJ). The densities , arc 



not the final stationary densities yet. After hitting the 
left boundary the shock reflects and changes its density 



to another value, u 2 



see Fig. [7] This reflected 



value 



wave, in its turn, hits the right boundary, reflects again 
and changes its density to u^ + ,v^ + . 

This process continues indefinitely and the sequence 
"} converges exponentially to the stationary 
„„it)iQat as n -> oo. The shock densities 
u n + i v n + f° r °dd n = 1,3,... (for even n = 2,4,...) be- 
long to the right reflection map defined by ur , vr (to the 
left reflection map defined by Ul,Vl), see 17|, [DJ for 
more details. In practice, it becomes harder and harder 
to observe reflections of high order since the respective 
shocks differ by an infinitesimal change of densities (see 



inset of FigJTJ)- The sequence {? 



"} , (unlike the 



fact of a presence of first order transition along the Phi 
path) is microscopic rates dependent, through the diffu- 
sion matrix B from (fTj). 

It is important to stress that precisely at the transition 
point the shock wave between the H — h state (matching 
the left boundary ul,vl) and the — h state with densi- 
ties , uj~ + is unbiased, meaning that there is a perfect 
balance between the respective currents: j u {uL,VL) = 
j u {ui + ,v~ + ), j v {u L ,v L ) = j v (u- + ,v~ + ). The fact that 
ci(ul,Vl) > (ci(% + , Vi + ) < 0) at the left (right) from 
shock discontinuity guarantees the stability of the unbi- 
ased shock G ++ /G |- . At one side of the phase transition 

this shock is biased to the right leading to G++ station- 
ary state, and at another side of the phase transition it is 

biased to the left leading to G |_ stationary state. This 

feature is essentially the same as in one-species systems, 
see [f| ■ Precisely at the phase transition line the unbiased 
shock performs a random walk between the boundaries. 
By averaging the local particle density over large times 
one samples configurations with the shock at all possible 
positions, which leads to a density profile with a linear 
slope, observed in Monte Carlo simulations (not shown 
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Table I: Sequential changes of the stationary state densities 
along a minimal path Phi ((5J from G++ to G |_. 



Table II: Sequential changes of the stationary state densities 
along a minimal path Ph2 from G++ to G_ + . 



for brevity). 

More complicated scenarios at first order transition 
points may be observed if we choose left and right bound- 
ary densities belonging to non-connected G domains (e.g. 
the left boundary belongs to G ++ and the right boundary 
belongs to G domain), and impose flat initial condi- 
tions matching one of the boundaries, e.g. the left bound- 
ary of G ++ type. Close to the phase transition G |_ — > 

G__ shown at FigfTU] along a single big Phi path, we 
shall see appearance of two shocks at the right boundary 

of + + / — h and of — h / type. For a while in the 

system there are two moving consecutive shocks. 

Stability of such a multi-shock is possible due to ex- 
istence of two conserved quantities, see 0. The first 
shock reaches the left boundary and reflects from it 
(with changed densities). Now we have two shocks which 
counter-propagate and collide at some point, forming a 

single shock of — V / type. The future stationary 

state is then determined by the direction of motion of this 
shock: for positive (negative) shock velocity the result- 
ing stationary state will be of — h (of type). Again, 

we see that the first order phase transition is caused 
by the direction of the bias of a shock connecting two 
states. The schematic space-time evolution of the above 
described scenario is shown in Fig|S] Note, however, that 
if we follow an adiabatic path, the initial conditions as we 
have imposed will not appear (close to a transition point 

G |- — > G the initial state, to be quasi-stationary, 

must be either of G |_ or of G__ type), and consequently 

at any time we shall see at most one shock in the system. 

This scenario of the first order phase transitions de- 
scribed above for two particle species is straightforwardly 
generalizable to an arbitrary number of species K. A dis- 
continuous phase transition of the type p (see Sect lIII[) 
along the minimal path described in Sec lIVl is caused by 
a shock Gx/Gy between a Gx state on the left from 
discontinuity and a Gy state on the right from discon- 
tinuity. State Gx has one extra positive characteristic 
speed (c p > 0) with respect to the Gy state (c p < 0), see 

The signs of all remaining characteristic speeds (but 
not the characteristic speeds themselves) are the same at 
both sides of the shock. At the transition point, the shock 
Gx/Gy is unbiased (has zero velocity) and its stability is 
guaranteed by the fact that c p > (c p < 0) on the left (on 
the right) from the discontinuity. Such a shock is called a 
p-shock in the PDE theory of conservation laws Q • Zero 
shock velocity signalizes equality of particle currents of 
all species at both sides of discontinuity. Consequently, 
the stationary current is continuous across the transition. 

At one side of the transition, the Gx / Gy shock is bi- 



ased to the right and then the G^-type stationary state 
prevails. At the other side of the transition the shock 
is biased to the left, resulting in the Gy-type stationary 
state. Thus, the density changes discontinuously at the 
transition point while the current is continuous (usually 
it has a cusp) across the transition point. The location of 
the transition point itself can be indicated on the respec- 
tive Phi path only approximately, as being inside the 
dashed-like decorated segment of it, see Figflja)) and 
Table 1. 

If the left and right boundary densities belong to non- 
connected G domains (but the number of positive char- 
acteristic speeds at the left boundary is larger than those 
at the right boundary) and if the initial state of the sys- 
tem matches one of the boundaries, a stable multiple 
shock can be observed. The stationary state of the sys- 
tem is then decided by a sequence of reflections from the 
boundaries (governed by the diffusion matrix B in U)), 
and interactions in the bulk between shocks (governed by 
the current-density relations j q {u\ 1 U2 1 ...,uk))- 



VI. RAREFACTION WAVE MECHANISM 
UNDERLYING SECOND ORDER BDPT 

If a Gx/Gy shock is stable (see proceeding Section), 
the inverse shock Gy/Gx is unstable and gives rise to 
a rarefaction wave which is a self-similar solution of ([1]), 
depending only on ratio £ = {x — xo)/t where xq is a 
position of its center, and t > 0. Let us argue that in the 
long-time limit t — > oo the stationary bulk density u s t a t 
generated by a rarefaction wave, has zero characteristic 
speed c p (u stQt ) = 0. By u stat we denote a set of bulk 
stationary densities {uf at , u^ 0,1 , u s ^ at }. We search for 
a solution of ((T|) in the form u(x, t) — /i(£). Substituting 
in ([T]) , and denoting the derivative with respect to £ with 
a prime, we obtain 

-ih' + ±(Dj)h' = ~0{e), (8) 

where the matrix (Dj)(h(£)) is the Jacobian of the flux 
(Dj) pq = dj p /du q . The above equation can be rewritten 

as 

(Dl)h' = th'+?Q. (9) 

In the limit t — > oo the 0(e)/t term vanishes, £ = 
(x — xo)/t — > for any finite x, and the above equa- 
tion reduces to (Dj^.^h' ~ 0, e.g. the solution is 
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an eigenvector of the flux Jacobian _Dj with zero eigen- 
value. Consequently, the (-Dj)t^oo = {Dj){u stat ) is a 
matrix with zero eigenvalue, i.e. u sta4 belongs to a sub- 
region Gyox with zero characteristic speed, situated "in 
between" Gy-type and Gx-type states. Such a subre- 
gion is the boundary between Gy-type and Gjf-type do- 
mains, a hyperplane of dimension K — 1 characterized 
by c p — 0. The respective rarefaction wave is called p- 
rarefaction wave Q,[8|. 

Arguments presented above and in Sec fVl imply a num- 
ber of consequences for the locations of continuous and 
discontinuous BDPTs, discussed below. 

Note that the scenario of a rarefaction wave governing 
long-time evolution may take place only if initial states 
Gy on the left (and Gx on the right) are supported by 
respective boundaries, meaning that left (right) bound- 
ary density is of Gy -type (of Gx type). Such a setting 
appears along a Ph2 path, see Sec lIVl and never appears 
along a Phi path. Inspecting a Ph2 path one finds that 
such a setting appears in the intermediate part of the Ph2 
path marked by dashed line in Fig|TJ starting as soon as 
the left boundary density crosses the c p = hyperplane 
(during step L) and finishing when the right boundary 
density crosses the c p — hyperplane (during step R). 
All along this intermediate Ph2 segment, the rarefaction 
wave governs the stationary state which stays "pinned" 
to the c p — hyperplane. Initial and final points of the 
segment are points where the pinning and depinning from 
the c p = hyperplane take place (see also Table 2). This 
conclusion is fully supported by numerical simulations. 
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Figure 9: Stationary densities pA, Pb, , pc , along a Phi path 

across domains G+++ — > G , versus running coordinate s 

along the path, represented by right boundary density of the 
first specie A along the path) for a three-chain model K = 3, 
in the torus setting where each chain have two chains- neigh- 
bours, see [l(|. A particle hops to the right neighbouring 
site with rate r n = 1 — where < n < 4 is the num- 
ber of particles on the adjacent chains, neighbouring to the 
departure and to the target sites. Parameters: 7 = 0.5. Ini- 
tial and final points of the path are INI= (0.2, 0.3, 0.4) and 
FIN = (1, 1, 1). Three discontinuous transitions, between the 
states G+++ — > G-++ — > G y — > G are clearly seen. 

Analogously, a shock wave leading to the discontin- 



uous phase transition is stable only if the left and the 
right boundary are of Gx- and of Gy-type (T5]) respec- 
tively. Such a setting always appears along a Phi path 
(the segment marked by bold dotted line in Figfl]). How- 
ever, for an existence of a stable unbiased shock, other 
conditions must be fulfilled, namely: (i) perfect balance 
between particle currents at both sides of discontinuity 
(ii) shock densities at both sides of discontinuity must 
form stable boundary layers with respective boundaries, 
i.e. to belong to respective reflection maps of Ul and ur 

m : 

Since the latter maps depend on the microscopic de- 
tails of the dynamic, see [13], [HI, we cannot locate pre- 
cisely the phase transition point, but deduce that it must 
be inside the unbiased shock-wave favourable segment 
marked by bold dotted line in FigfT] 

On the other hand, the unbiased shock-wave favourable 
setting never appears along a Ph2 path. Therefore, dis- 
continuous changes in stationary densities described by 
our shock wave scenario, cannot happen along Ph2 path. 
Consequently, any Ph2 path p s tat(s) in physical region 
is always continuous, see Fig. Q] (b). Reciprocally, a 
favourable setting for stable rarefaction wave formation 
never appears along a Phi path. Therefore, a state with 
c p = 0, governed by a stable rarefaction wave, cannot be 
observed along any Phi path. Consequently, since ini- 
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Figure 10: Location of stationary densities along the single 
Phi path (filled circles) and single Ph2 path (open triangles) 
from G++ to G — domain, for two-chain model with 7 = 0.5. 
Evolution direction is marked by arrows. Crosses show the 
initial and final points. 

tial and final stationary states pjni and pfin belong to 
different regions with c p > and c p < 0, at least one 
discontinuous change must happen along any Phi path, 
see Fig. Q](b). 

It should be clear from our reasoning that one can con- 
struct other, more complicated paths in parameter space, 
along which one can observe the same phenomenon of 
discontinuous or continuous phase transitions. Any path 
in parameter space connecting points pini and pfin in 
different G-regions, and not containing segments favour- 
ing rarefaction waves, will result in discontinuous phase 
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transitions in physical region (i.e. will be Phi-like). Re- 
ciprocally, any path not containing segments favouring 
shock-waves, (and containing therefore favorable bound- 
ary settings for stable rarefaction waves formation, will 
show only continuous phase transitions (i.e. will be Ph2 
-like). Further examples are given below. 



VII. SPECIAL PATHS FOR SEQUENCES OF 
BDPTS 

Special sequences of phase transitions in system with 
K species can be observed along rather simple paths. 

i) Phi and Ph2 paths between disjoint G\ and 
Gy domains. 

Along any single (not composite) Phi-like path ([5]) 
connecting arbitrary disjoint G regions, a sequence of 
first order phase transitions will be observed, provided 
that the initial state has more positive characteristic 
speeds than the final state. The latter condition makes 
existence of stable shocks (governing first order tran- 
sitions) possible and leads to observation of as many 
discontinuous transitions, as the number of hyperplanes 
c p = separating the initial and the final state. The 
existence of such a path was pointed out in Q. For ex- 
ample, to observe all qualitatively different first order 
transitions, we can take a single Phi path from initial 
point with all positive characteristics /9, n , G G+..+ to a fi- 
nal point with all negative characteristics p final £ 
Along such a path, K first order transitions will be ob- 
served, see Figs llOI9l and the path marked by squares in 
Fig lllf a). In particular, Fig|S]corresponds to multi-chain 
model with K = 3 and shows respectively three discon- 
tinuous transitions in stationary density along the Phi 
path. The model has product stationary states which 
allows to compute the particle fluxes and consequently 
characteristic velocities analytically as functions of par- 
ticle densities. 

By computing the characteristic velocities along the 
path in physical region p s tat{s) for K = 3 we find that 
across each discontinuous transition just one character- 
istic velocity changes sign, this confirming the first-order 
phase transition scenario described in Sec. IVl 

Similarly, a single Ph2 path from initial to final state 
which belong to disjoint G regions allows to observe 
the sequence of all continuous transitions between these 
states. As examples, see Figs. [9l [10] (see also the path 
marked by squares in Fig lTTT b) of the next section for a 
more complex model) . Note that it is important that that 
the initial state has more positive characteristic speeds 
than the final state, 
ii) Fully matching path. 

Another special path is a path where left and right 
boundary densities are equal all along from the initial 
till the end path point. It is clear that in this case we 
will not observe any phase transitions because there will 
be always a perfect match of the bulk density with the 
boundaries. Consequently, this path in parameter space 




(b) 

Figure 11: Stationary densities of right and left-movers 
(u s tat, Vstat) = Pstat for bidirectional traffic model from Monte 
Carlo simulations for a system with 300 sites, along various 
Phi paths (Panel (a)) and Ph2 paths (Panel (b)). Refer- 
ence Initial, Middle and Final points are marked by crosses 
INI, MID and FIN). Lines where one characteristic veloc- 
ity is zero (Go+,G-o) are obtained numerically. Evolution of 
the pstat(s) along a path is marked by arrows, dotted arrows 
mark discontinuous transitions. Few data points outside the 
arrow-marked paths result from finite size effects. The sym- 
metry of the Figure with respect to the line y = x is due to the 
left-right symmetry of the model and the points INI, FIN . 
Parameters: h = 0.5. Panel (a): Squares mark pstat along 
a Phi path which goes directly from the initial to the final 
point INI— s> FIN. Triangles and circles mark pstat along two 
consecutive Phi paths INI-> MID, MID^ FIN. Initial, 
Middle and Final boundary rates, corresponding to points 
INI, MID, FIN are (a = 1 - p3 = 0.1, A = 1 - B = 0.78); 
(a = l-f3 = A = l-B = 0.9) and (A = 1 - B = 0.1, a = 
1 — ft = 0.78). Along all paths, the boundary rates a, /3,A,B 
are changed by linear interpolation law. Panel (b): The 
same as Panel (a), for respective Ph2 paths. Note that the 
densities along the direct Ph2 path (squares) go through the 
weak hyperbolic point ci = C2 = 0, marked by M, see also 
discussion at the end of Sec lVIII Initial, Middle and Final 
boundary rates: (a = 1 - p = 0.1, A = 1 - B = 0.75), 
(a = l-P = A = l-B = 0.95), and (A = 1 - B = 0.1, a = 
1 — f3 = 0.75) respectively. 
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must contain all triple points where the hyperplanes of 
second order and first order transitions merge together. 
For K = 1 such a triple point is a point {p* , p* ) where the 
characteristic speed vanishes j'(p*) = 0. The nature and 
topology of the parameter space in vicinity of these triple 
points (lines, hypersurfaces) will be discussed elsewhere. 



VIII. MODEL FOR A BIDIRECTIONAL 
TRAFFIC ON A NARROW ROAD 

In the previous sections we concentrated our attention 
on solvable models with analytic flux functions and strict 
hypcrbolicity. For generic (not intcgrablc) models, how- 
ever, the analytic flux function is typically unknown, as 
well as exact relation between boundary rates and ef- 
fective reservoir boundary densities. In the following we 
show how even in this case it is still possible to construct 
Phi-like or Ph2-like paths along which a given BDPT 
type (discontinuous or continuous) can be observed. As 
an example, we consider the case of a two-way traffic 
model on a narrow road (see bottom panel of Fig. 3). 

Models of bidirectional traffic have been widely stud- 
ied in the literature and appear in several contexts see 
e.g. [llj]. Our system consists of two chains, containing 
particles hopping in opposite direction: a particle hops in 
preferred direction with constant rate 1 and hard core ex- 
clusion like in TASEP, but is slowing down when it meets 
an upcoming particle (an obstacle) in front on the adja- 
cent lane: in this case the rate of hopping is exp(— h), 
where a positive constant h measures the interlane in- 
teraction, see Figj3] A similar model, but with periodic 
boundary conditions, was considered in [l9|. We choose 
the boundary rates as follows: if the target site is vacant, 
a right-moving particle can enter with rate a (ae~ h ) if 
the adjacent to the target site is empty (is occupied by 
an upcoming particle). At the other end, a particle can 
leave with rate j3. For the left moving particles, the en- 
trance and exit rates are respectively A(Ae~ h ) and B. 
Note, that the model has the left-right symmetry. Since 
the model is not solvable, the analytical expression for 
the flux j(u,v,h) is not known for any nonzero h. Nei- 
ther we know the exact relation between the boundary 
rates and the effective boundary densities. 

Phi- and Ph2-like paths, however, can be constructed 
straightforwardly. From the physical meaning of the 
characteristic velocities (e.g. velocities with which small 
perturbations of the homogeneous state propagate) [l3| 
we conclude that a stationary state with small density 
of right moving particles or right moving holes realized 
e.g. for a = 1 — < 1 and B = 1 — A < 1, has all 
positive characteristic velocities and therefore it must be 
in the G++ region. By left-right symmetry, a stationary 
state with small density of left moving particles or holes 
will belong to the G region ( the respective bound- 
ary rates are attainable from G++ rates by exchanging 
a A, f3 <==> B). Finally, a state with small density of 
right movers on one lane and small density of left movers 



on another lane, realized by a, A, 1 — j3, 1 — B <C 1 or 
l-a,l-4,/3,B<Cl, belongs the G y region. Proceed- 
ing along Phi (Ph2) paths in parameter space between 

regions G ++ —> G |_, and G |_ — > G__, one expects 

to see the occurrence of first (second) order phase transi- 
tions as described above. This is precisely what we obtain 
from Monte Carlo simulations of the two-way model, see 
FigE] 

We can also build direct Phi and Ph2 paths between 

G++ — > G as described in Sec lVIII i, see square data 

points in Fig|TT] Note that along a direct Ph2 path (see 
FigHHb)), the pinning/depinning of stationary densi- 
ties to the line with c p = occurs only once, due to 
a presence of a special point (or region) M in the mid- 
dle where the lines c\ = and C2 = intersect. Such 
a point where two characteristic velocities coincide (the 
so- called weakly hyperbolic point), makes possible a con- 
tinuous passage from G++ to G domain. It is worth 

to note that according to the numerical study of Jiang 
et al. [l9[ restricted to the case of periodic boundary 
conditions and equal particle densities, the steady state 
current along the symmetric line u sta t = Vstat develops a 
plateau, leading in periodic system to phase separation. 
Such a non-analiticity in the stationary current suggests 
that the region M in the middle of Figfrijb) is a segment 
rather that a single point. 

It is quite remarkable that even in this, rather special 
situation with non-analytic current-density dependence, 
our predictions about discontinuity/continuity of phase 
transitions along Phl/Ph2 paths remain robust. We also 
remark that the presence of the region M can be neglected 
as long as our paths are situated far enough from it. The 
systematic study of an influence of a weakly hyperbolic 
point on BDPTs will be done elsewhere. 



IX. CONCLUSIONS 

In this paper we have classified the basic phase transi- 
tions which can be observed in multi-species driven sys- 
tems with open boundaries. We have shown that the 
splitting of the physical region into domains with differ- 
ent signs of characteristic speeds, and hyper-surfaces sep- 
arating these regions where one of characteristic speeds 
vanishes, plays a fundamental role in this classification. 
Adiabatic paths in the parameter space, defined by the 
particle densities of each specie at the left and right 
boundary reservoirs, along which we surely observe dis- 
continuous or continuous transition of a desired type, or 
a desired sequence of BDPTs, have been explicitly con- 
structed. The details of the microscopic dynamics and 
the geometry of the models are not important for our 
qualitative BDPTs scenarios to occur, as far as several 
conditions listed at the beginning of Sec|TI] are fulfilled. 
We expect therefore our results to be valid for a broad 
class of particle models with several interacting particle 
species. In particular, our examples were systems of par- 
ticles obeying hard-core exclusion rule, but this is not 
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required as far as some interaction making the flux func- 
tion nonlinear will be present. 

Mathematically, our study has been focused mainly to 
models with analytic flux function and strict hyperbol- 
icity i.e. models with Jacobian matrices which have dis- 
tinct eigenvalues in all the physical region. An example 
of a weakly hyperbolic model with non-analytic flux and 
phase separation, however, was considered in Sec I Villi It 
is remarkable that even for this model the general validity 
of our approach has been confirmed. An interesting prob- 
lem for the future would be to test the predictions of our 
analysis on more complicated models, like those showing 
symmetry breaking, hysteresis and ergodicity breaking 
phenomena. 
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